%% task6_2_2_4.m

clear;
[y, Fs] = audioread("..\wav\fmt.wav");
% Fs = 1;
% 去除前后静音片段
Mth = 1e-3;         % -60dB
fmt_begin = find(y>Mth, 1);
fmt_end = find(y<Mth, 1, 'last');
fmt = y(fmt_begin:fmt_end);

% audiowrite("..\wav\fmt_new.wav", fmt, Fs);
figure
t = 0 : 1/Fs : (length(fmt)-1)/Fs;
plot(t, fmt)
hold on

% 绘制包络
% 自带函数，并不好用
% [up,lo] = envelope(fmt, 262, "peak");
% plot(t,up,t,lo,'linewidth',1.5)

% ampd尝试
% loc = ampd(up);
% save("..\data\fmt_loc.mat", "loc");
% load("..\data\fmt_loc.mat", "loc");
% xline((loc-1)/Fs)

% 另一种方法，切隔音频
[pks, locs] = enve(fmt, 4);
[up, lup] = findpeaks(pks);
[lo, llo] = findpeaks(-pks);

Data = [locs(lup), up;locs(llo), -lo];
Data = sortrows(Data);
temp = [];
j = 1;
for i = 1 : length(Data)-1
    if Data(i+1,2)-Data(i,2) > 0.04
        temp(j) = Data(i+1,1);
        j = j+1;
    end
end
sou = [temp(1)];
j = 2;
for i = 2 : length(temp)
    if temp(i) - temp(i-1) > 1088
        sou(j) = temp(i);
        j = j+1;
    end
end
sou = [sou(1:9),sou(11:length(sou))];       % 手动删除第10个点，吉他拨弦太慢，目前没有合适算法
plot((locs-1)/Fs, pks);
xline(sou/Fs, 'c');
title('时域图像')
% xline(locs(lup)/Fs, 'm');
% xline(locs(llo)/Fs, 'c');

% fft_plot(fmt(sou(1):sou(2)), 'Fig1');

% 通过频谱图计算每个音的频率
[S,F,T,P] = spectrogram(fmt, 2048, 2000, 4000, Fs);
figure(2);
surf(T,F,10*log10(P),'edgecolor','none');
view(0,90);
title('时频分析');
xlabel('t/s');
ylabel('f/Hz');
axis([0 16 0 4000]);

% 重新排序
[sortP,I] = sort(P,'descend');
sortf = F(I);
% 频率检测
Pth1 = 1e-4;     % 1st threshold of power
Pth2 = 1e-5;     % 2nd threshold of power

funds = cell(size(sortf,2),1);  % to save fundamental results

for i = 1:size(sortf,2)
    % fundamental = [];
    col = sortP(:,i);
    if max(col)>Pth1
        f = sortf(1:5,i);       % 5 most powerful components
    else
        f = [];
    end
    if isempty(f)
        f = sortf((col>Pth2));
        if isempty(f)
            f = [];
        end
    end
   
    fundamental = ffilter(f);   % detect fundamental component
    % fundamental may equals [1312,656], as f is power-sorted
    % so we must sort fundamental and ffilter again
    fundamental = ffilter(sort(fundamental,'ascend'));
    
    funds{i} = fundamental;
end

% 匹配到最近的乐音，以小字组a为基准，220Hz
ns = cell(length(T),1);
for i = 1:length(T)
    ns{i} = round(12*log2(funds{i}/220));
end

% 去除噪音以及非线性谐波
% 注意到时间向量的间隔T(2)-T(1)=0.006, 假设人耳的分辨率为resolution,
% 那么稳定的音应当持续resolution/0.006个时间间隔以上, 按此条件滤除时间较短的噪音:
resolution = 0.15;
combo = round(resolution/(T(2)-T(1)));
for i = 2:length(T)-combo
    if ~isequal(ns{i},ns{i-1}) && ~isequal(ns{i},ns{i+combo})
        ns{i} = ns{i-1};
    end
end

% 解析每个音的时值
nt = cell(1,2);
nt{1,1} = ns{1}; nt{1,2} = T(1);

for i = 2:length(T)
    if isequal(ns{i},nt{end,1}) || isequal(ns{i},[])
        nt{end,2} = T(i);
    else
        nt = [nt;cell(1,2)];
        nt{end,1} = ns{i}; nt{end,2} = T(i);
    end
end

for i = 1:size(nt,1)-1
    nt{end-i+1,2} = nt{end-i+1,2} - nt{end-i,2};
end

[song,BPM,key]=interpreter(nt);

disp(song);

save("..\data\task6_2_2_4.mat", "sortf", "sortP", "T", "Fs", "fmt");


% 迭代计算包络
function [pks, locs] = enve(data, iter)
    if iter > 1
        [pks, l1] = findpeaks(data);
        [pks, l2] = enve(pks, iter-1);
        locs = l1(l2);
    else
        [pks, locs] = findpeaks(data);
    end
end

function fft_plot(data, name)
    L = length(data);
    Y = fft(data);
    % 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
    P2 = abs(Y/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    
    f = 8000*(0:(L/2))/L;
    figure
    plot(f,P1) 
    title(strcat(name, " 频谱分析"));
    xlabel("f (Hz)");
    ylabel("|P1(f)|");
end

function fundamental = ffilter(f)
% 在输入的一系列频率值中, 找出可能的基频
    fundamental = [];
    for j = 1:length(f)
        if isempty(fundamental)
            fundamental = [fundamental,f(j)];
        else
            ratio = f(j)./fundamental;
            inrange = (ratio<(round(ratio)*1.05)) + (ratio>(round(ratio)*0.95));    % 近似整数倍
            if all(inrange~=2)
                fundamental = [fundamental,f(j)];
            end
        end
    end
end

function [song,BPM,key] = interpreter(nt)
%[song,BPM,key] = interpreter(nt)
%输入:
%   <n-by-2 cell> nt: 第一列为相对于220Hz差几个半音, 第二列为当前音的长度
%输出:
%   <n-by-2 cell> song: C大调简谱, 第二列为节拍数
%   <float> BPM: 每分钟节拍数
%   <char> key: 'C'
    key = 'C';
    t = cell2mat(nt(:,2));
    tpb = min(t);
    BPM = round(60/tpb);
    song = cell(size(nt,1),2);
    for i = 1:size(nt,1)
        nnote = [];
        n = nt{i,1};
        for j = 1:length(n)
            nnote = [nnote,floor(n(j)/12)*7+cnotetable(mod(n(j),12))];
        end
        song{i,1} = nnote;
        song{i,2} = round(nt{i,2}/tpb);
    end
end

function nnote = cnotetable(n)
    switch n
        case 0
            nnote = -1;
        case 1
            nnote = -0.5;
        case 2
            nnote = 0;
        case 3
            nnote = 1;
        case 4
            nnote = 1.5;
        case 5
            nnote = 2;
        case 6
            nnote = 2.5;
        case 7
            nnote = 3;
        case 8
            nnote = 4;
        case 9
            nnote = 4.5;
        case 10
            nnote = 5;
        case 11
            nnote = 5.5;
    end
end
